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Abstract. We apply a Monte-Carlo method to the two flavor Polyakov loop extended 
Nambu and Jona-Lasinio (PNJL) model. In such a way we can go beyond mean field 
calculations introducing fluctuations of the fields. We study the impact of fluctuations on the 
thermodynamics of the model. We calculate the second derivatives of the thermodynamic grand 
canonical partition function with respect to the chemical potential and present a comparison 
with lattice data also for flavor non-diagonal susceptibilities. 



1. Introduction 

Results of QCD thermodynamics from lattice computations are reproduced surprisingly well 
with a quasiparticle model, an extension of the Nambu - Jona-Lasinio model with inclusion of 
Polyakov loop dynamics (the PNJL model) at the mean field level (D0E1HE1E1E1I8]. A better 
understanding of these results requests the investigation of fluctuations in the PNJL model. 
This can be done by numerical simulations of the thermodynamics using standard Monte-Carlo 
techniques (MC-PNJL). The advantage over other methods is that Monte-Carlo calculations 
automatically incorporate fluctuations to all orders. 

We present new results of such computations and discuss the role and importance of 
fluctuations beyond mean field. It turns out that the temperature dependence of the second 
non-diagonal derivative of the thermodynamic grand canonical partition function with respect 
to quark chemical potentials is particularly sensitive to such fluctuations. While the mean field 
PNJL predicts a vanishing coefficient, our calculation agrees well with available lattice data. 

2. The PNJL partition function 

The Euclidean action of the two-flavor PNJL model, at finite baryon and isospin chemical 
potential, is given by (9j [7] 





where V is the Nf = 2 doublet quark field and ttiq = diag(m u , is the quark mass matrix. 
The quark chemical potential matrix is defined as fl = diag(/x + fj,i,fj> — fj,i), where [±i is the 
isospin isovector chemical potential. In the PNJL model quarks interact with a background 
color gauge field A 4 = iA , where A = 5^gA^t a with A% G SU(3) C and t a = X a /2. The field 
A^ is related to the traced Polyakov loop by: 

$ = T^-tr c L with L = exp U j£ drAA and /3 = — . (2) 
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In the Polyakov gauge, the matrix L is given in a diagonal representation 

L = exp(i((j) 3 X 3 + </> 8 A 8 )). (3) 
The dimensionless effective fields 03 and 0§ introduced here are identified with the Euclidean 

(3) (8) 

gauge fields in temporal direction divided by the temperature, A\ jT and A\ /T. These two 
fields are a parametrisation of the diagonal elements of SU(3) C . 

The effective potential U models the confinement-deconfinement transition in pure gauge 
QCD on mean-field Ginzburg-Landau level. In this paper we consider an ansatz for the effective 
potential given in [51 [10] motivated by the SU(3) Haar measure: 

U(<S>,$*,T) = _ 1 Q ^ $ * $ + b ^ ln ^ _ 6$ * $ + 4 ^*3 + _ 3 ($*$)2] ? ( 4 ) 
where the temperature dependent prefactors are given by 

a(T) = a + ai(^)+a 2 (^) 2 and b(T) = b 3 ^ . (5) 

The parameters are chosen such that the critical temperature of the first order transition is 
equal to Tq = 270 MeV, as given by lattice calculations, and that <!>*, — > 1 as T — > oo. 

After performing a bosonization of the PNJL Lagrangian introducing the auxiliary fields a 
and 7r the partition function can be rewritten as 



Z = N J Vcj)VaViTexp[-S[a(x),7r(x),<j)(x)}} 

= M j ' VcpVaVve^-TrHS- 1 ]- [3 J d 3 x(l4(<f>,fi) + (^t?L))]. (6) 

with the inverse quark propagator 

s -\ _ ( + {n + m - iA 4 )7o - M \ 
V *T57r -@ + (fi- M- iA 4 )>y -M J ' 

and M = tuq — a. 

In order to introduce fluctuations in the PNJL model we consider the total volume of our 
system as composed of a certain number of subregions, V = J2i V%, and require that the fields 
are completely correlated if they lie in the same subvolume, and decorrelated otherwise |11] . 
In this way the fields can fluctuate in the sense that they can have different values from one 
subvolume to another. Physically this means that our fields are only weakly dependent on the 
position and can be considered constant inside subvolumes of a given dimension (Fig. [TJ. 

From this assumption it follows that the action S[a(x), 7r(cc), <j>{x)\ can be decomposed into 
a sum of contributions related to the i-th volume 
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Figure 1. Approximation of the space dependence of the fields: the behavior of the field (blue 
line) is approximated in our model using a mean value (red lines) on subvolumes. In this way 
we permit small fluctuations of the field. 



S[a(x),ir(x),(f>{x)} = ^ S[a(x), 7r(x), <l>{x)]\ SeV . = ^ vr, </>]. (7) 
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The last equality emphasizes the strong correlation of the fields inside each subvolume. We 
assume therefore that the fields are constant for all x £ Vi, hence, the partition function becomes 

Z = ]J ( jV(j)iVaiV7riexp[-Si(a,TT,(f))]). (8) 
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Computing the partition function Q we introduce fluctuations in the PNJL model. 
3. Monte-Carlo method 

In order to evaluate Q, we adopt standard Monte-Carlo techniques. Using the Metropolis 
algorithm we start from a random configuration C and suggest a new configuration C that will 
be accepted with probability 

p = mm{l,exp{-S[C'])/exp{-S[C})}. (9) 

Writing Eq. Q, we point out that the fields involved in each subvolume are independent 
from the fields in all the other volumes. As a consequence, if we update a given configuration 
moving only the fields in the k-th. subregion the action changes according to 

n 

s(0i,...,0' fc) ...&O = + 

n 

S((j>i,...,(f)k,..;(l>n) = ^2S((f>i) 

i=l 

s(<j)i, ...,(/)' k , ...,<j) n ) - s((pi, ...,<p k , ■■ -An) = s(4>' k ) - s((j>k) (io) 

this means that the Metropolis algorithm acts independently on each single volume. If we thus 
assume that every Vi will be of a fixed size, Vi = V, Eq. QSJ) can be rewritten in momentum 
space as 

Z = (J V<j>VaVirexp[l3V^J2 J ^3 Tl MS _1 (^n,p)] 

2 2 

-pv(u(<P,(3) + (^^-))]) n n-oo. (11) 

The calculation of the partition function using Monte-Carlo can be performed by generating a 
set of configurations for a given subvolume. We only need to fix the size of this volume. Here 
we adopt the standard temperature dependence of the volume used in lattice calculations, 

1 /V 3 

a= KT - V = N - a3 =N^- <12) 



where a is the lattice spacing, Nt is the number of lattice sites in the Euclidean time direction 
and N s the number of lattice sites in the space dimensions. It follows that 

V oc k/T 3 , (13) 

where k will be conveniently chosen. In this paper we perform calculations at different values 
of k in order to understand the impact of k on thermodynamic quantities. In the next sections 
we will see how this Monte-Carlo version of the PNJL model (MC-PNJL) works both in pure 
gauge and in the Nj = 2 case. 

4. Pure Gauge MC-PNJL 

In this section we apply the Monte-Carlo approach to pure gluodynamics, i.e. to the Polyakov 
loop sector of the model. Contributions due to fluctuations in this case are very small, although 
we prefer to start with this simple case, because it is instructive for the understanding of the 
method. 

We start with a generic configuration that in our system corresponds to A% and ^4s at a 
given temperature. Using Monte-Carlo we generate points in the configuration space close to 
the minimum of the action. In particular applying cooling we will reach exactly the minimum. 
In Fig. [2] we can see the difference between cooling and Metropolis in the case of T = T c = 0.27 
GeV. For cooling (left picture) the action must decrease, and in few steps we fall into one of 
the three minima of the potential. On the other hand, using the full Metropolis the action can 
also slightly increase. As a consequence, the space of the available configurations is much larger 
(right picture). 




Figure 2. Comparison between the cooling (left) and Metropolis (right) procedure applied to 
the Polyakov loop effective potential at T = T c = 0.27 GeV. 

To fix the parameters of the effective potential we have computed the pressure, energy and 
entropy density at different temperatures including fluctuations, trying to fit the available lattice 
data (Fig. [3] left). In the right panel of Fig. [3] we compare the Polyakov loop obtained using 
the mean-field method (dashed line) with our results (red points) and the lattice data (black 
points) we can see that there is practically no difference between the two methods: both agree 
well with lattice data. 

Once established how this method works in the case of the pure gluodynamics we can move 
to the more interesting case of the two flavors MC-PNJL model. 




Figure 3. Left: fit of Lattice data [12] for pressure, energy and entropy density. Right: 
Comparison of the Polyakov loop dynamics with the data in |13j 



5. Two flavors MC-PNJL model 

The starting point studying the thermodynamics for Nf = 2 is the partition function ( |11[ ). The 
degrees of freedom in this case are the ^3 and A$ components of the gauge field, and the two 
bosonic fields a and ir. As explained before, evaluating the path inte gral , we need to fix the size 
of the volume V as a function of the temperature. According to Eq. ( 13 ) we only have to fix the 
value of k. Since the value of this parameter is not known a priori, in this paper we consider 
four different values, k = 95,135,500,5000. In this way we can study the dependence of the 
observables as a function of the volume size. 



5.1. Chiral and deconfinement transitions 

Since our aim is the investigation of the QCD phase diagram, the first thing we need to study is 
how the chiral and deconfinement transitions are affected by the introduction of fluctuations in 
the model. This can be achieved evaluating the chiral condensate and the Polyakov loop for the 
different volumes and comparing that with the mean-field result. This comparison is presented 
m Fig. |4j the presence of fluctuations does not modify the behavior of the Polyakov loop, the 
five different sets of data overlap perfectly. For the chiral condensate we can notice only small 
dependence on the temperature below the critical temperature. In analogy to the pure gauge 
case we can conclude for Nf = 2 that the order parameters of the transitions are not affected 
by the presence of fluctuation. 

This last observation is indeed important when we move to the case of finite chemical 
potential, because this means that we can use mean-field approximation to determine the critical 
temperature: corrections due to the presence of fluctuations are negligible, and this permits to 
overcome the sign problem which is not present in mean field calculations. 



6. Taylor expansion at finite chemical potential 

The Taylor expansion approach to QCD thermodynamics at non-zero quark chemical potential 
represents one of the possibilities to overcome the sign problem on the lattice: instead of 
performing explicit calculations at \i q 7^ one expands the thermodynamic potential in a Taylor 
series around zero chemical potential, 

n=0 



(14) 
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Figure 4. Chiral condensate and Polyakov loop: dependence on the volume size. The red lines 
are the mean field calculation, while moving from the blue to the green line we go from k = 95 
to k = 5000. 
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1 d n lnZ 



(15) 



m=o n\VT s d(n/T) n m=o 

and n even due to the CP symmetry of the system. 

The comparison between lattice data and MC-PNJL calculations for these coefficients 
represents a crucial test of the model. In particular the second non-diagonal coefficient is the 
most important because, as explained in the introduction, to obtain a non-vanishing result for 
we must take into account fluctuations of the fields. We know that the size of fluctuations 
in our model depends on the size of our correlation volume, therefore, as before, we evaluate the 
Taylor coefficients for different values of the parameter k. 

To perform the calculations using Monte-Carlo, we need to represent the derivative of the 
logarithm of the partition function in an explicit form. To do that, we start from the definition 
of our partition function, 

-V 



Z[T,ii,ix I ]= VfVAexp - In det M[T, M , ^i, f, A] - S g [ A] 



(16) 



From this expression we have 
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Z[T, (i, m] T 



J VfVA 



d In det M[T,n,fj, I} f, A] 



dfj, 



-S[T,n,m,f,A] 



V i d In det M[T,/x, m,f, A] 



(17) 



Proceeding in the same way for the second derivative we obtain for the C2 and isovector 
coefficients 

1 T 2 d 2 \nZ[T^,^J,A] 



2 VT 3 



Q..2 



1 / d 2 In det M[T, /x, fjj, /, A] 

tF 2 \ 



(18) 



where we have neglected the expectation values which involves odd derivatives, due to the 
symmetry property of the determinant. 

The definition of the flavor diagonal and non-diagonal coefficients we are interested in is given 

by 

uu _ c n + c n ud _ c n ~ C„ 

n ~ 4 n ~ 4 

In Fig. [5] we show the results for df 1 and c% d . For the diagonal coefficient we can see that the 
different choices of k affect the matching with the lattice data (black points) only slightly. The 
picture changes totally when looking at the non-diagonal coefficient: in this case, at increasing 
volume the signal reduces. This is exactly what we expect because we know that with growing V 
the fluctuations responsible for the non vanishing expectation value disappear. Thus we recover 
the vanishing mean-field result (dashed light- blue line). 



(19) 




Figure 5. Second order diagonal and non-diagonal coefficients cJf 4 and cV; . Black points are 
lattice data from p3|; the colored lines show the model predictions (k increase from blue to red) 



It is quite surprising, though, that for the smallest values of k we can reproduce the lattice 
data with reasonable accuracy. This seems to indicate that the PNJL model can capture indeed 
the dynamics of the fluctuations at the origin of the non- vanishing c| . 



7. Conclusions 

In this work we applied standard Monte-Carlo techniques in order to go beyond mean-field in 
the Polyakov loop extended Nambu - Jona-Lasinio model calculations. Fluctuations introduced 
in this way strongly depend on a new parameter that fixes the minimal volume V as defined in 
Sec. [2j Inside this volume the fields are totally correlated while they are totally de-correlated on 
distances larger than V 1 ^ 3 . Studying the QCD thermodynamics we found that the introduction 
of fluctuations does not alter the behavior of the traced Polyakov loop in the pure gauge sector 
nor the behavior of the Polyakov loop and chiral condensate in Nf = 2. On the other hand these 
fluctuations are crucial in the evaluation of the Taylor expansion coefficients of the pressure. In 
particular, using the MC-PNJL we can reproduce very well the lattice data for the second 
non-diagonal coefficient which vanishes in mean field approximation. 

Collecting all this results we have that, the MC-PNJL model on the one hand includes cor- 
rectly the dynamics associated with confinement and chiral symmetry breaking, on the other 
hand we can reproduce quite well lattice data. So it runs for as a good model in the investigation 
of the full phase diagram. 
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